Elranatamab & Teclistamab IPTW Analysis

1 SETUP

1.1 Open libraries

library(gtsummary)
library(dplyr)
library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(readxl)
library(janitor)
library(ggplot2)
library(tidycmprsk)
library(ggsurvfit)
library(tidymodels) 
library(sjPlot)
library(patchwork)
library(bstfun)
library(ComplexUpset)

library(writexl)
library(rms)

library(timeROC)
library(rms)
library(pec)
library(prodlim)
library(riskRegression)

library(plotly)
library(ggtext)

library(WeightIt)
library(cobalt)
library(survey)

2 DATA WRANGLING

2.1 Read in MM datasets

db_Kansas = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 - due 5-1-25 (de-identified).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  filter(!is.na(age_at_bs_ab_initiation)) %>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Kansas, "Kansas_elra_data.xlsx")

db_CC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/CC bsAb Consortium Submission CCF 6-11-25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_CC, "CC_elra_data.xlsx")

db_HUMC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/BsAb Consortium Data HUMC.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_HUMC, "HUMC_elra_data.xlsx")

db_MCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data MCC-deidentified.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_Huntsman = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/Utah_Huntsman_bsAb Consortium Data FINAL_5.4.25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc)/1000,
    baseline_alc = as.double(baseline_alc)/1000,
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Huntsman, "Huntsman_elra_data.xlsx")


db_FHCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/FHCC Bispecifics (5.5.2025).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    date_of_infection_80 = ymd(date_of_infection_80),
    baseline_anc = as.double(baseline_anc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_UTSW = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 1.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3")) %>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles)
  )

# Write to Excel
#write_xlsx(db_UTSW, "UTSW_elra_data.xlsx")


db_MDACC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/MDACC BITE ASH 2025.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_MDACC, "MDACC_elra_data.xlsx")


db_Roswell = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 5.8.2025_NN.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_Roswell, "Roswell_elra_data.xlsx")

# Combine all into a single dataframe
combined_db <- bind_rows(
  list(
    HUMC = db_HUMC,
    MCC = db_MCC,
    Huntsman = db_Huntsman,
    FHCC = db_FHCC,
    UTSW = db_UTSW,
    MDACC = db_MDACC,
    Roswell = db_Roswell,
    CC = db_CC,
    Kansas = db_Kansas
  ),
  .id = "site"
) %>%
  mutate(
    prior_bcma_directed_therapy_yes_no = ifelse(!is.na(most_recent_bcma_agent_name),1,0),
    bs_ab_name = case_when(
      bs_ab_name == 1 ~ "tec",
      bs_ab_name == 3 ~ "elra",
      .default = NULL
    ),
    max_crs_grade_0_5 = ifelse(is.na(max_crs_grade_0_5), 0,max_crs_grade_0_5), ## Set NAs to 0
    max_icans_grade = ifelse(is.na(max_icans_grade), 0, max_icans_grade ) ## Set NAs to 0
  )

# Convert all datetime columns to Date
combined_db_clean <- combined_db %>%
  mutate(across(where(lubridate::is.POSIXt), as.Date))

# Write to Excel
write_xlsx(combined_db_clean, "elra_tec_data.xlsx")


# Write to Excel
# write_xlsx(combined_db_clean %>% filter(prior_bcma_directed_therapy_yes_no != prior_bcma_directed_therapy), "elra_data_discordant.xlsx")

## Several institutions do not have any cases of elranatamab

# db_UABMC = read_excel(path = "Data/bsAb Consortium Data 5.8.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Stanford = read_excel(path = "Data/De-identified Stanford bsAb Consortium Data Template 4.30.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_LCI = read_excel(path = "Data/bsAb Consortium Data LCI.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_COH = read_excel(path = "Data/bsAb Consortium Data Template 2.11.2025_COHv2.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_MUSC = read_excel(path = "Data/USMM bsAb Consortium Data Template 2.11.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Mayo = read_excel(path = "Data/Bispecific_consortium_20250521 NoPHI.xlsx", sheet = "bsAb_DATA") %>%
#   filter(`bsAb name` == 3)%>%
#   clean_names()

2.2 Calculate CAR-HEMATOTOX score and risk level

combined_db <- combined_db %>% 
  mutate(
    # baseline_anc = as.numeric(baseline_anc),
    # baseline_hgb,
    # baseline_pl_ts,
    # baseline_crp_prior_to_bs_ab_initiation,
    # baseline_ferritin_prior_to_bs_ab_initiation,
    HT_Score = 
      (baseline_anc <= 1.2) +
      (baseline_hgb <= 9.0) +
      (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) +
      (baseline_crp_prior_to_bs_ab_initiation >= 3.0) +
      (baseline_ferritin_prior_to_bs_ab_initiation >= 650 & baseline_ferritin_prior_to_bs_ab_initiation <= 2000) +
      (baseline_pl_ts <= 75) + ## Add an extra point for Plt <=75
      (baseline_ferritin_prior_to_bs_ab_initiation > 2000), ## Add an extra point for ferritin >2000
    
    # Classify the risk based on the score
    HT_Risk = factor(
      case_when(
        HT_Score >= 2 ~ "HThigh",
        baseline_pl_ts <=75 ~ "HThigh", ## Equivalent to at least 2 points
        baseline_ferritin_prior_to_bs_ab_initiation > 2000 ~ "HThigh", ## Equivalent to at least 2 points
        (baseline_anc <= 1.2) + (baseline_hgb <= 9.0) + (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) >=2 ~ "HThigh",
        HT_Score <2 ~ "HTlow",
        .default = "NA"
      ),
      levels = c("HTlow", "HThigh")
    )
)
write_xlsx(combined_db %>% select(site, patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate, baseline_anc, baseline_hgb, baseline_pl_ts, baseline_crp_prior_to_bs_ab_initiation, baseline_ferritin_prior_to_bs_ab_initiation, HT_Score, HT_Risk), "HT_score.xlsx")

2.3 Identify severe infections

combined_db <- combined_db %>% 
  mutate(
    severe_infection = ifelse(infection_severity_73 %in% 1 | infection_severity_76 %in% 1 | infection_severity_79 %in% 1| infection_severity_82 %in% 1,1,0),
    first_severe_infection_date = case_when(
      infection_severity_73 %in% 1 ~ date_of_infection_74,
      infection_severity_76 %in% 1 ~ date_of_infection_77,
      infection_severity_79 %in% 1 ~ date_of_infection_80,
      infection_severity_82 %in% 1 ~ date_of_infection_83,
      .default = NULL
    ),
    first_severe_infection_type = case_when(
      infection_severity_73 %in% 1 ~ number_infection_1_bacterial_viral_fungal,
      infection_severity_76 %in% 1 ~ number_infection_2_bacterial_viral_fungal,
      infection_severity_79 %in% 1 ~ number_infection_3_bacterial_viral_fungal,
      infection_severity_82 %in% 1 ~ number_infection_4_bacterial_viral_fungal,
      .default = NULL
    )
    
    
  )

3 PROPENSITY SCORES

3.1 Regression analysis for treatment allocation

theme_gtsummary_compact()

variables <- c("elra", "Age", "Female sex", "ECOG 2+", "EMD", "Hgb", "LDH", "Prior BCMA")

data1 <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  transmute(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    elra = ifelse(bs_ab_name == "elra", 1, 0),
    "ECOG 2+" = ifelse(ecog_ps_at_bs_ab_initiation_0_5 >=2,1,0),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

uv_tab <- tbl_uvregression(
    data1[c(variables)],
    method = glm,
    y = "elra",
    method.args = list(family = binomial),
    exponentiate = TRUE
  ) %>%
  sort_p()


mv_tab<-glm(elra ~ Age + `Female sex` + `ECOG 2+` + EMD + Hgb + LDH + `Prior BCMA`, data=data1, family=binomial) %>%
  tbl_regression(exponentiate=TRUE)

tbl_merge(list(uv_tab, mv_tab), tab_spanner = c("**Univariate**", "**Multivariable**"))

Characteristic	Univariate				Multivariable
	N	OR1	95% CI1	p-value	OR1	95% CI1	p-value
Hgb	371	1.12	1.01, 1.25	0.029	1.14	1.02, 1.28	0.026
EMD	371	1.69	0.95, 2.96	0.069	1.71	0.96, 3.04	0.068
Age	371	1.02	1.00, 1.04	0.087	1.01	0.99, 1.04	0.3
Female sex	371	1.43	0.93, 2.21	0.11	1.43	0.92, 2.24	0.11
Prior BCMA	371	0.76	0.49, 1.17	0.2	0.78	0.50, 1.23	0.3
ECOG 2+	371	1.19	0.75, 1.88	0.4	1.27	0.78, 2.07	0.3
LDH	371	1.00	1.00, 1.00	0.8	1.00	1.00, 1.00	>0.9
1 OR = Odds Ratio, CI = Confidence Interval

3.2 IPTW

IPTW.data <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  mutate(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    bsAb = bs_ab_name,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 %in% c("0","1") ~"0-1",
      ecog_ps_at_bs_ab_initiation_0_5 ==2 ~"2",
      ecog_ps_at_bs_ab_initiation_0_5 >=3 ~"3+",
      .default = NULL),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

## Use the WeightIt function to generate propensity scores and weights
W.out <- weightit(bsAb ~ Age + `ECOG` + EMD + Hgb + LDH + `Prior BCMA`,
                  data=IPTW.data, 
                  method = "glm",    ## Methods: glm, gbm, cbps, npcbps, ebal, optweight, super, bart, energy
                  estimand = "ATE",
                  stabilize = TRUE)

bal.tab(W.out, stats = c("m", "v"), thresholds = c(m=0.25), binary = "std")

## Balance Measures
##                Type Diff.Adj     M.Threshold V.Ratio.Adj
## prop.score Distance   0.0050 Balanced, <0.25      0.9733
## Age         Contin.   0.0090 Balanced, <0.25      0.9595
## ECOG_0-1     Binary  -0.0121 Balanced, <0.25           .
## ECOG_2       Binary   0.0152 Balanced, <0.25           .
## ECOG_3+      Binary  -0.0027 Balanced, <0.25           .
## EMD          Binary  -0.0034 Balanced, <0.25           .
## Hgb         Contin.  -0.0186 Balanced, <0.25      1.0974
## LDH         Contin.  -0.0276 Balanced, <0.25      1.0961
## Prior BCMA   Binary  -0.0119 Balanced, <0.25           .
## 
## Balance tally for mean differences
##                     count
## Balanced, <0.25         9
## Not Balanced, >0.25     0
## 
## Variable with the greatest mean difference
##  Variable Diff.Adj     M.Threshold
##       LDH  -0.0276 Balanced, <0.25
## 
## Effective sample sizes
##              elra    tec
## Unadjusted 123.   248.  
## Adjusted   113.38 242.39

summary(W.out)

##                   Summary of weights
## 
## - Weight ranges:
## 
##         Min                                  Max
## elra 0.4745 |---------------------------| 2.3134
## tec  0.7811      |-------------|          1.7418
## 
## - Units with the 5 most extreme weights by group:
##                                         
##          123    102     81     91      9
##  elra 1.5458 1.5896  1.653 1.7607 2.3134
##          169     71    274    258    271
##   tec 1.4257 1.4453 1.5215 1.6783 1.7418
## 
## - Weight statistics:
## 
##      Coef of Var   MAD Entropy # Zeros
## elra       0.292 0.222   0.040       0
## tec        0.152 0.115   0.011       0
## 
## - Effective Sample Sizes:
## 
##              elra    tec
## Unweighted 123.   248.  
## Weighted   113.38 242.39

love_plot <- love.plot(W.out,
          binary = "std",
          line = TRUE, 
          abs = TRUE,
          thresholds = c(m=0.2),
          sample.names = c("Unweighted", "IPTW"),
          title = "A)",
          var.order = "unadjusted",
            var.names = c(
              prop.score = "Distance"
          ))+ theme(plot.title = element_text(hjust = 0))
love_plot

weight_df = data.frame(W.out$weights, W.out$treat)

IPTW.data$weights <- W.out$weights
  
IPTW.data$ps <- W.out$ps


write.csv(data1, "IPTW.dataset.csv")

3.3 Scatterplot of weights

Weights_plot <- ggplot(
  IPTW.data, aes(x=ps, y = weights, color = bsAb)) + 
  geom_point(shape = 16) +
  theme_classic() + 
  scale_y_continuous(breaks = seq(0,4 , by = 0.5)) +
  labs(x = "Propensity score", y = "Weights") + 
  ggtitle("B)")
Weights_plot

4 TABLES

4.1 Site contributions

theme_gtsummary_compact()

combined_db %>%
  transmute(
    "Site" = site,
    bs_ab_name,
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    )
  ) %>%
  bold_labels()

Characteristic	elra N = 1301	tec N = 2771
Site
MCC	53 (41%)	72 (26%)
MDACC	12 (9.2%)	57 (21%)
CC	20 (15%)	47 (17%)
UTSW	10 (7.7%)	40 (14%)
FHCC	7 (5.4%)	24 (8.7%)
Huntsman	9 (6.9%)	22 (7.9%)
Kansas	4 (3.1%)	13 (4.7%)
HUMC	10 (7.7%)	1 (0.4%)
Roswell	5 (3.8%)	1 (0.4%)
1 n (%)

4.2 Unweighted and IPTW baseline characteristics

theme_gtsummary_compact()

data1 = IPTW.data %>%
  transmute(
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    bs_ab_name,
    weights
  )

table_unweighted <- data1 %>% select(-weights) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels()

table_IPTW <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  ) %>%
  tbl_svysummary(
    by = bs_ab_name,
    missing = "no",
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"),
    include = c("Age (years)", "Female sex","Race","ECOG","LDH (U/L)","Hgb (g/dL)","Baseline ferritin (ng/mL)", "Heavy chain","True EMD","≥1 HRCA","Prior LOTs","Prior ASCT","Prior CAR-T","Triple refractory","Penta refractory","Trial eligible","Prior BCMA-directed therapy", "Received full dose")
  ) %>%
  add_p() %>%
    bold_labels()
  
tbl_merge(list(table_IPTW, table_unweighted), tab_spanner = c("**IPTW**", "**Unweighted**"))

Characteristic	IPTW			Unweighted
	elra N = 1231	tec N = 2481	p-value2	elra N = 1233	tec N = 2483	p-value4
Age (years)	70 (64, 76)	69 (62, 76)	0.9	71 (39 to 95)	68 (35 to 88)	0.060
Female sex	66 (54%)	114 (46%)	0.2	68 (55%)	115 (46%)	0.11
Race			0.2			0.15
White	94 (79%)	187 (76%)		95 (81%)	187 (76%)
Black	17 (14%)	50 (20%)		16 (14%)	51 (21%)
AAPI	3 (2.5%)	6 (2.3%)		3 (2.5%)	5 (2.0%)
American Indian	2 (1.8%)	0 (0%)		2 (1.7%)	0 (0%)
Other	2 (2.0%)	4 (1.8%)		2 (1.7%)	4 (1.6%)
ECOG			>0.9			0.5
0	23 (18%)	49 (20%)		21 (17%)	50 (20%)
1	62 (50%)	119 (48%)		59 (48%)	121 (49%)
2	28 (23%)	59 (24%)		29 (24%)	60 (24%)
3+	10 (8.2%)	20 (8.2%)		14 (11%)	17 (6.9%)
LDH (U/L)	223 (178, 280)	208 (167, 290)	0.4	222 (175, 277)	208 (167, 291)	0.5
Hgb (g/dL)	9.90 (8.50, 11.40)	9.80 (8.40, 11.50)	0.9	10.00 (8.70, 11.70)	9.60 (8.25, 11.20)	0.027
Baseline ferritin (ng/mL)	346 (140, 1,250)	379 (125, 1,219)	>0.9	293 (137, 832)	393 (128, 1,254)	0.3
Heavy chain			0.001			<0.001
IgG	72 (58%)	109 (44%)		72 (59%)	107 (43%)
IgA	26 (21%)	43 (17%)		26 (21%)	42 (17%)
IgM	1 (0.7%)	1 (0.4%)		1 (0.8%)	1 (0.4%)
IgD	2 (1.5%)	0 (0%)		2 (1.6%)	0 (0%)
None	23 (19%)	95 (38%)		22 (18%)	98 (40%)
True EMD	20 (16%)	39 (16%)	>0.9	26 (21%)	34 (14%)	0.067
≥1 HRCA	81 (70%)	151 (65%)	0.3	81 (70%)	152 (65%)	0.4
Prior LOTs	6.00 (5.00, 8.00)	6.00 (4.00, 8.00)	0.4	6.00 (4.00, 8.00)	6.00 (4.00, 8.00)	>0.9
Prior ASCT	66 (54%)	141 (57%)	0.6	61 (50%)	144 (58%)	0.12
Prior CAR-T	53 (43%)	90 (36%)	0.2	50 (41%)	94 (38%)	0.6
Triple refractory	112 (91%)	224 (91%)	0.9	112 (91%)	226 (91%)	>0.9
Penta refractory	61 (49%)	123 (49%)	>0.9	61 (50%)	126 (51%)	0.8
Trial eligible	24 (19%)	55 (22%)	0.5	26 (21%)	53 (21%)	>0.9
Prior BCMA-directed therapy	67 (54%)	133 (54%)	>0.9	60 (49%)	138 (56%)	0.2
Received full dose	114 (93%)	241 (97%)	0.059	114 (93%)	241 (97%)	0.045
1 Median (Q1, Q3); n (%)
2 Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment
3 Median (Min to Max); n (%); Median (Q1, Q3)
4 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

5 STACKED BAR GRAPH

5.1 FIGURE 1A: Toxicity, IPTW

data1 = IPTW.data %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    CRS.any = ifelse( max_crs_grade_0_5 != 0, 1, 0),
    ICANS.any = ifelse( max_icans_grade != 0,1,0),
    
    CRS.two = ifelse( max_crs_grade_0_5 == 2 | max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.two = ifelse( max_icans_grade ==2 | max_icans_grade == 3 | max_icans_grade == 4,1,0),
    
    CRS.three = ifelse( max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.three = ifelse( max_icans_grade == 3 | max_icans_grade == 4,1,0),
    weights
  )


weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )


combined_data <- rbind(
  IPTW.data %>%
    count(bs_ab_name, max_crs_grade_0_5, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      Category = "CRS"
    ),
  IPTW.data %>%
    count(bs_ab_name, max_icans_grade, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )


# Any grade toxicity
CRS.pval <- as.numeric(svychisq(~CRS.any+BsAb, weighted_table)$p.value)
ICANS.pval <- as.numeric(svychisq(~ICANS.any+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G2+ toxicity
CRS.two.pval <- as.numeric(svychisq(~CRS.two+BsAb, weighted_table)$p.value)
ICANS.two.pval <- as.numeric(svychisq(~ICANS.two+BsAb, weighted_table)$p.value)

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G3+ toxicity
CRS.three.pval <- as.numeric(svychisq(~CRS.three+BsAb, weighted_table)$p.value)
ICANS.three.pval <- as.numeric(svychisq(~ICANS.three+BsAb, weighted_table)$p.value)

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )




toxplot_IPTW <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(palette = "Pastel1") +
  #guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 115)) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("A) IPTW")

toxplot_IPTW

5.2 FIGURE 1B: Toxicity, unweighted

combined_data <- rbind(
  combined_db %>%
    count(bs_ab_name, max_crs_grade_0_5) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "CRS"
    ),
  combined_db %>%
    count(bs_ab_name, max_icans_grade) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
    
  )




## Any-grade toxicity
CRS.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 !=0,1,0) ,
  test = combined_db$bs_ab_name
)
CRS.pval <- chisq.test(CRS.d.test$results,CRS.d.test$test, correct = FALSE)$p.value

ICANS.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade !=0,1,0) ,
  test = combined_db$bs_ab_name
)
ICANS.pval <- chisq.test(ICANS.d.test$results,ICANS.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )


## G2+ toxicity
CRS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 == 2 | combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.two.pval <- chisq.test(CRS.two.d.test$results,CRS.two.d.test$test, correct = FALSE)$p.value

ICANS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 2 | combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.two.pval <- chisq.test(ICANS.two.d.test$results,ICANS.two.d.test$test, correct = FALSE)$p.value

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )

## G3+ toxicity
CRS.three.d.test <- data.frame(
  results = ifelse( combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.three.pval <- fisher.test(CRS.three.d.test$results,CRS.three.d.test$test)$p.value

ICANS.three.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.three.pval <- chisq.test(ICANS.three.d.test$results,ICANS.three.d.test$test, correct = FALSE)$p.value

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )



toxplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 100)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = 1) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("B) Unweighted")

toxplot_unweighted

5.3 ***FIGURE 1: IPTW and unweighted toicity

FIGURE1 <- toxplot_IPTW / toxplot_unweighted

ggsave("svg/FIGURE1.svg", FIGURE1, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE1.tiff", FIGURE1, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE1.jpeg", FIGURE1, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE1.jpeg")

5.4 FIGURE 2A: Response, IPTW

data1 = IPTW.data %>%
  filter(!is.na(best_response)) %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR",1,0),
    VGPR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR",1,0),
    CR = ifelse(best_response == "sCR" | best_response == "CR",1,0),
    weights
  )

weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )



data1 <-   IPTW.data %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response, wt=weights) %>%
  group_by(bs_ab_name) %>%
  mutate(
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR", "SD","PD")),
    bs_ab_name = as.factor(bs_ab_name),
    Category = "Response rate"
  )



Response.pval <- as.numeric(svychisq(~ORR+BsAb, weighted_table)$p.value)
CR.pval <- as.numeric(svychisq(~CR+BsAb, weighted_table)$p.value)
VGPR.pval <- as.numeric(svychisq(~VGPR+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

responseplot_IPTW <- ggplot(data1) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR")) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR"| best_response == "PR"),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(
    palette = "Pastel1",
    direction = -1,
    breaks = rev(levels(data1$best_response))
  ) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels  ) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("A) IPTW")

responseplot_IPTW

5.5 FIGURE 2B: Response, unweighted

data1 <- combined_db %>%
  filter(!is.na(best_response)) %>%
  transmute(
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" |best_response == "PR",1,0),
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    Category = "Response rate"
  )

combined_data <- combined_db %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response) %>%
  group_by(bs_ab_name) %>%
  transmute(
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    n,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR"]),
    Category = "Response rate"
  )

## ORR
Response.d.test <- data.frame(
  results = ifelse(data1$ORR,1,0),
  test = data1$bs_ab_name
)
Response.pval <- chisq.test(Response.d.test$results,Response.d.test$test, correct = FALSE)$p.value

## VGPR
VGPR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR" | data1$best_response == "VGPR",1,0),
  test = data1$bs_ab_name
)
VGPR.pval <- chisq.test(VGPR.d.test$results,VGPR.d.test$test, correct = FALSE)$p.value

## CR
CR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR",1,0),
  test = data1$bs_ab_name
)
CR.pval <- chisq.test(CR.d.test$results,CR.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )


responseplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR")) +
  geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
             data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"),
             position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = -1) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("B) Unweighted")

responseplot_unweighted

5.6 ***FIGURE 2: IPTW and unweighted response

FIGURE2 <- responseplot_IPTW / responseplot_unweighted

ggsave("svg/FIGURE2.svg", FIGURE2, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE2.tiff", FIGURE2, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE2.jpeg", FIGURE2, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE2.jpeg")

6 SURVIVAL ANALYSIS

6.1 Length of follow-up (Reverse KM, months)

data1 <- IPTW.data %>%
  transmute(
    bs_ab_name,
    Death = ifelse(death_yes_no,1,0),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Reverse_death = ifelse(death_yes_no == 1, 0,1),
    weights
  )

reverse_km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1)
reverse_km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     356 356    215   13.4    12.5    14.9

reverse_km_OS_elra <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "elra"))
reverse_km_OS_elra

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "elra"), weights = weights)
## 
##    6 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     117 117   74.5   7.47    6.03       9

reverse_km_OS_tec <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "tec"))
reverse_km_OS_tec

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "tec"), weights = weights)
## 
##    9 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     239 239    140   17.1    14.9    20.1

6.2 Stratified by BsAb

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     117 117   42.8     NA     8.6      NA
## bs_ab_name=tec      239 239   98.6   21.5    16.4      NA

med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    12 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     118 118   68.1   3.77    2.57    7.17
## bs_ab_name=tec      241 241  151.1   8.50    7.10   11.30

med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    3 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      62  60.2   20.6   11.7    6.13      NA
## bs_ab_name=tec      118 120.7   56.1   13.7    9.27    22.3

med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.3 ***FIGURE 3: PFS/OS/DOR, stratified by trial elibility

FIGURE3 <- wrap_plots(
  (OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
  (PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
  (DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
  ncol = 1  # Force vertical stacking
)

ggsave("svg/FIGURE3.svg", FIGURE3, device = "svg",
       width = 6.5,
       height = 8.5,
       units = c("in"))

ggsave("jpeg/FIGURE3.jpeg", FIGURE3, device = "jpeg",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE3.tiff", FIGURE3, device = "tiff",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE3.jpeg")

6.3.1 Schoenfeld residuals and PH test

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    elra_ind = ifelse(bs_ab_name == "elra", 1, 0),  # 1 = elra, 0 = tec
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

## ---- Cox models ----
cox_OS <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ elra_ind, weights = weights, data = data1)

cox_PFS <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ elra_ind, weights = weights, data = data1)

## ---- PH diagnostics ----
zph_OS  <- cox.zph(cox_OS)
zph_PFS <- cox.zph(cox_PFS)

# Global PH p-values (same as elra_ind right now)
p_OS  <- signif(zph_OS$table["GLOBAL",  "p"], 3)
p_PFS <- signif(zph_PFS$table["GLOBAL", "p"], 3)

par(mfrow = c(1, 2))

plot(
  zph_OS,
  main = paste0(
    "Schoenfeld residuals – OS\n",
    "Global PH test: p = ", p_OS
  )
)

plot(
  zph_PFS,
  main = paste0(
    "Schoenfeld residuals – PFS\n",
    "Global PH test: p = ", p_PFS
  )
)

par(mfrow = c(1, 1))

6.4 High risk subsets

6.4.1 True EMD

data1 <- IPTW.data %>%
  filter(EMD == 1) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(EMD == 1) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    4 observations deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      24 18.2   4.86     NA    5.23      NA
## bs_ab_name=tec       32 37.4  21.05   14.1   10.43      NA

med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    1 observation deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      25 18.8   8.27   6.27    2.57      NA
## bs_ab_name=tec       34 39.5  27.45   7.93    4.67    12.5

med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    1 observation deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      11  8.76   2.39     NA    5.27      NA
## bs_ab_name=tec       15 18.05   8.69   6.93    5.63      NA

med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.4.2 ECOG 2+

data1 <- IPTW.data %>%
  filter(ecog_ps_at_bs_ab_initiation_0_5 >= 2) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(ecog_ps_at_bs_ab_initiation_0_5 >= 2) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    4 observations deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      42 37.9   18.8   4.87    2.27      NA
## bs_ab_name=tec       74 75.8   45.8   9.57    5.27    18.8

med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    2 observations deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      43 38.6   27.9   1.63    1.07   10.27
## bs_ab_name=tec       75 76.7   56.7   4.43    2.10    8.73

med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      15 12.9   6.02   11.7     1.0      NA
## bs_ab_name=tec       30 31.3  16.87   13.3     8.3      NA

med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.4.3 Hgb <10 g/dL

data1 <- IPTW.data %>%
  filter(baseline_hgb <10) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(baseline_hgb <10) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    8 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      55  63.3   31.7   4.87    3.17      NA
## bs_ab_name=tec      134 125.3   71.5  12.83    9.57    18.6

med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    4 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      56  63.9   43.7   1.63    1.10    6.27
## bs_ab_name=tec      137 128.3   97.8   3.70    2.17    7.40

med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    1 observation deleted due to missingness 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      26 28.6   13.6   6.17     4.2      NA
## bs_ab_name=tec       51 47.5   27.4   9.90     7.0      NA

med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.4.4 Age ≥80 years

data1 <- IPTW.data %>%
  filter(age_at_bs_ab_initiation >=80) %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  filter(age_at_bs_ab_initiation >=80) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS

## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      18 15.0   11.1    5.2    3.47      NA
## bs_ab_name=tec       30 32.9   10.6   17.7    9.53      NA

med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##                 records    n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      18 15.0   11.1   3.73    1.63      NA
## bs_ab_name=tec       30 32.9   14.8   8.63    6.07      NA

med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR

## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra       9  8.13   5.47   4.20    1.00      NA
## bs_ab_name=tec       19 20.84   7.68   6.63    5.07      NA

med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.5 ***FIGURE 4: PFS/OS/DOR, stratified by Age ≥80 yrs

FIGURE4 <- wrap_plots(
  (OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
  (PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
  (DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
  ncol = 1  # Force vertical stacking
)

ggsave("svg/FIGURE4.svg", FIGURE4, device = "svg",
       width = 6.5,
       height = 8.5,
       units = c("in"))

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).

ggsave("jpeg/FIGURE4.jpeg", FIGURE4, device = "jpeg",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"))

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).

ggsave("tiff/FIGURE4.tiff", FIGURE4, device = "tiff",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).

knitr::include_graphics("jpeg/FIGURE4.jpeg")

7 REGRESSION

7.1 ORR, ≥CR, ≥VGPR, PFS, OS, DOR, IFS, CRS G1+, ICANS G1+, steroid use (IPTW)

options(gtsummary.add_p_footnote = FALSE, gtsummary.add_estimate_to_reference_rows = FALSE)
theme_gtsummary_compact()


preds0 <- c("Response", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds1 <- c("CR/sCR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds2 <- c("BsAb", "Baseline Hgb",  "Baseline LDH", "ECOG 2+","True EMD", "Prior LOTs", "Prior BCMA")
preds3 <- c("ICANS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds4 <- c("CRS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds5 <- c("Steroid use", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

preds6 <- c("VGPR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

data0 <- IPTW.data %>%
  transmute(
    weights,
    BsAb = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    HT_Risk,
    best_response,
    "Infection" = ifelse(infection == 1, 1,0),
    "CrCl <40" = ifelse(cr_cl_prior_to_step_up_dose_number_1 <40, 1,0),
    "CR/sCR" = ifelse(best_response %in% c("CR","sCR"),1,0),
    "VGPR" = ifelse(best_response %in% c("CR","sCR", "VGPR"),1,0),
    "Response" = ifelse(best_response %in% c("CR","sCR", "VGPR", "PR"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Steroid use" = ifelse(steroid_use_as_tx==1,1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG 2+" = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "ICANS 1+" = case_when(
      is.na(max_icans_grade) ~ NA,
      max_icans_grade %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "CRS 1+" = case_when(
      is.na(max_crs_grade_0_5) ~ NA,
      max_crs_grade_0_5 %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    prior_bcma_directed_therapy_yes_no,
    Days.since.last.BCMA = case_when(
      is.na(prior_bcma_directed_therapy_yes_no) ~ NA,
      date_of_last_exposure_to_prior_bcma_agent > date_of_bs_ab_step_up_d1 ~ NA,
      prior_bcma_directed_therapy_yes_no == 1 ~ as.numeric(as.Date(date_of_bs_ab_step_up_d1) - as.Date(date_of_last_exposure_to_prior_bcma_agent))
    ),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    "Baseline Hgb" = baseline_hgb,
    "Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/100,
    "Baseline ferritin" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/1000,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
      ),
    date_of_bs_ab_step_up_d1,
    bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes,
    days_to_infection_or_last_contact = as.numeric(
      pmin(as.Date(date_of_infection_74), as.Date(date_of_last_contact), na.rm = TRUE) - as.Date(date_of_bs_ab_step_up_d1)
    ),
    infection_or_death = ifelse(!is.na(date_of_infection_74) | death_yes_no == 1,1,0),
    IVIg = ifelse(ivig_use_yes_1_no_0 == 1, 1, 0),
    ALPS = ifelse(baseline_hgb >10 & baseline_ldh_prior_to_bs_ab_initiation <300,1,0),
    baseline_pl_ts,
    baseline_anc,
    baseline_hgb,
    baseline_crp_prior_to_bs_ab_initiation,
    baseline_ferritin_prior_to_bs_ab_initiation,
    CAR.HEMATOTOX = baseline_pl_ts * -0.05 + baseline_anc * -0.03 + baseline_hgb * -0.20 + as.numeric(baseline_crp_prior_to_bs_ab_initiation) * 0.15 + as.numeric(baseline_ferritin_prior_to_bs_ab_initiation) * 0.001
    ) %>%
  mutate(
    "BCMA exposure status" = factor(case_when(
      #prior_bcma_directed_therapy_yes_no == 0 ~ "  Never",
      Days.since.last.BCMA < 365 ~ "  <1 yr",
      Days.since.last.BCMA >= 365 ~ "  ≥1 yr"
    ), c("  ≥1 yr","  <1 yr")),
    `CAR-HEMATOTOX risk` = factor(
      case_when(
        CAR.HEMATOTOX < 0 ~ "Low",
        CAR.HEMATOTOX <= 5 ~ "Intermediate",
        CAR.HEMATOTOX > 5 ~ "High"
      ),
      levels = c("Low", "Intermediate", "High")
  )
  )

data1 <- data0 %>% filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) ## Exclude patients who received bsAb as bridging

data2 <- data1 %>% 
  filter(best_response %in% c("PR","VGPR","sCR", "CR")) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)


## Response
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds0)),
  method = glm,
  y = `Response`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `Response` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds0), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**") 

table_response <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## VGPR or better
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds6)),
  method = glm,
  y = `VGPR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `VGPR` ~ BsAb +`True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds6), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_VGPR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## sCR/CR
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds1)),
  method = glm,
  y = `CR/sCR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `CR/sCR` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds1), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_sCR.CR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## OS

uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.DLC, data1$Death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.DLC, Death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data1,  weights = weights) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_os <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


## PFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_pfs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## DOR
uv_tab <- tbl_uvregression(
  data2[c(preds2)],
  method = coxph,
  y = Surv(data2$Days.to.death.or.relapse.or.DLC, data2$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data2$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data2) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_dor <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## IFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$days_to_infection_or_last_contact, data1$infection_or_death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
 ) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  days_to_infection_or_last_contact, infection_or_death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ifs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## ICANS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds3)),
  method = glm,
  y = `ICANS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `ICANS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds3), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ICANS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## CRS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds4)),
  method = glm,
  y = `CRS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `CRS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds4), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_CRS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## Steroid use
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds5)),
  method = glm,
  y = `Steroid use`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-glm( `Steroid use` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds5), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_steroid <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


# Combine all tables into one
final_table <- tbl_stack(
  list(
    table_response %>% modify_header() %>% modify_footnote(everything() ~ NA, abbreviation = TRUE) %>%
      modify_footnote(estimate_1 ~ "Estimate = Odds ratio or hazard ratio", abbreviation = TRUE),
    table_VGPR,
    table_sCR.CR,
    table_os,
    table_pfs,
    table_dor,
    table_ifs,
    table_ICANS,
    table_CRS,
    table_steroid
  ),
  group_header = c(
    "Response",
    "≥VGPR",
    "CR/sCR",
    "Overall Survival",
    "Progression-free Survival",
    "Duration of Response",
    "Infection-free Survival",
    "ICANS 1+",
    "CRS 1+",
    "Steroid use"
  )
)
  
final_table

Characteristic	Univariate				Multivariable
	N	Estimate1	95% CI	p-value	Estimate	95% CI	p-value
Response
BsAb	370
Tec		—	—		—	—
Elra		0.86	0.55, 1.35	0.5	0.86	0.52, 1.42	0.6
Baseline Hgb	370	1.46	1.29, 1.65	<0.001	1.34	1.18, 1.53	<0.001
Baseline LDH	370	0.86	0.77, 0.94	0.002	0.86	0.77, 0.95	0.004
ECOG 2+	370	0.39	0.25, 0.61	<0.001	0.42	0.26, 0.69	<0.001
True EMD	370	0.91	0.52, 1.63	0.8	0.94	0.50, 1.78	0.9
Prior LOTs	370	0.87	0.79, 0.94	0.001	0.91	0.82, 1.01	0.087
Prior BCMA	370	0.56	0.36, 0.85	0.008	0.50	0.29, 0.87	0.014
≥VGPR
BsAb	370
Tec		—	—		—	—
Elra		0.62	0.40, 0.95	0.031	0.57	0.35, 0.93	0.026
Baseline Hgb	370	1.38	1.23, 1.55	<0.001	1.32	1.17, 1.49	<0.001
Baseline LDH	370	0.84	0.73, 0.93	0.003	0.84	0.73, 0.93	0.004
ECOG 2+	370	0.52	0.33, 0.81	0.004	0.57	0.34, 0.93	0.026
True EMD	370	0.63	0.35, 1.10	0.11	0.56	0.30, 1.04	0.071
Prior LOTs	370	0.88	0.80, 0.95	0.003	0.95	0.86, 1.06	0.4
Prior BCMA	370	0.47	0.31, 0.72	<0.001	0.41	0.24, 0.70	0.001
CR/sCR
BsAb	370
Tec		—	—		—	—
Elra		0.70	0.44, 1.10	0.12	0.68	0.40, 1.12	0.13
Baseline Hgb	370	1.44	1.29, 1.63	<0.001	1.39	1.23, 1.58	<0.001
Baseline LDH	370	0.86	0.74, 0.96	0.019	0.87	0.75, 0.98	0.044
ECOG 2+	370	0.51	0.31, 0.81	0.006	0.61	0.36, 1.03	0.068
True EMD	370	0.66	0.35, 1.20	0.2	0.62	0.31, 1.19	0.2
Prior LOTs	370	0.89	0.81, 0.97	0.014	0.95	0.85, 1.06	0.4
Prior BCMA	370	0.56	0.37, 0.86	0.009	0.52	0.31, 0.89	0.018
Overall Survival
BsAb	352
Tec		—	—		—	—
Elra		1.41	0.96, 2.07	0.083	1.46	0.97, 2.20	0.073
Baseline Hgb	352	0.73	0.66, 0.81	<0.001	0.77	0.69, 0.86	<0.001
Baseline LDH	352	1.13	1.06, 1.20	<0.001	1.15	1.10, 1.20	<0.001
ECOG 2+	352	2.25	1.61, 3.16	<0.001	1.92	1.32, 2.78	<0.001
True EMD	352	1.34	0.88, 2.03	0.2	1.28	0.82, 2.00	0.3
Prior LOTs	352	1.06	0.99, 1.13	0.091	1.04	0.97, 1.11	0.3
Prior BCMA	352	1.09	0.77, 1.54	0.6	1.20	0.81, 1.78	0.4
Progression-free Survival
BsAb	355
Tec		—	—		—	—
Elra		1.43	1.05, 1.94	0.021	1.52	1.11, 2.09	0.009
Baseline Hgb	355	0.77	0.71, 0.83	<0.001	0.81	0.74, 0.89	<0.001
Baseline LDH	355	1.18	1.08, 1.28	<0.001	1.17	1.08, 1.27	<0.001
ECOG 2+	355	1.63	1.23, 2.16	<0.001	1.43	1.04, 1.95	0.027
True EMD	355	1.04	0.74, 1.47	0.8	0.94	0.67, 1.32	0.7
Prior LOTs	355	1.10	1.05, 1.16	<0.001	1.06	1.00, 1.12	0.056
Prior BCMA	355	1.58	1.18, 2.10	0.002	1.63	1.17, 2.27	0.004
Duration of Response
BsAb	177
Tec		—	—		—	—
Elra		1.49	0.87, 2.56	0.14	1.50	0.85, 2.63	0.2
Baseline Hgb	177	0.84	0.74, 0.95	0.007	0.83	0.72, 0.96	0.014
Baseline LDH	177	1.05	0.90, 1.23	0.5	1.04	0.89, 1.21	0.7
ECOG 2+	177	1.05	0.64, 1.72	0.8	0.89	0.52, 1.52	0.7
True EMD	177	1.48	0.75, 2.93	0.3	1.36	0.68, 2.72	0.4
Prior LOTs	177	1.01	0.89, 1.15	0.8	1.02	0.90, 1.17	0.7
Prior BCMA	177	1.08	0.67, 1.74	0.8	1.06	0.65, 1.75	0.8
Infection-free Survival
BsAb	358
Tec		—	—		—	—
Elra		1.25	0.93, 1.67	0.13	1.24	0.92, 1.67	0.2
Baseline Hgb	358	0.87	0.81, 0.94	<0.001	0.90	0.83, 0.97	0.010
Baseline LDH	358	1.10	1.05, 1.15	<0.001	1.09	1.05, 1.14	<0.001
ECOG 2+	358	1.70	1.30, 2.22	<0.001	1.52	1.14, 2.02	0.004
True EMD	358	1.16	0.83, 1.62	0.4	1.10	0.79, 1.54	0.6
Prior LOTs	358	0.98	0.93, 1.03	0.4	0.98	0.93, 1.04	0.5
Prior BCMA	358	0.85	0.65, 1.11	0.2	0.99	0.73, 1.35	>0.9
ICANS 1+
BsAb	370
Tec		—	—		—	—
Elra		0.87	0.49, 1.52	0.6	0.91	0.49, 1.65	0.8
Baseline Hgb	370	0.83	0.72, 0.95	0.008	0.84	0.72, 0.98	0.024
Baseline LDH	370	1.05	0.95, 1.15	0.3	1.04	0.93, 1.14	0.5
ECOG 2+	370	3.03	1.77, 5.21	<0.001	2.44	1.39, 4.31	0.002
True EMD	370	1.40	0.69, 2.69	0.3	1.31	0.62, 2.62	0.5
Prior LOTs	370	0.88	0.77, 0.98	0.034	0.88	0.76, 1.01	0.069
Prior BCMA	370	0.61	0.36, 1.03	0.067	1.01	0.53, 1.92	>0.9
CRS 1+
BsAb	370
Tec		—	—		—	—
Elra		0.53	0.34, 0.82	0.004	0.55	0.35, 0.86	0.010
Baseline Hgb	370	0.94	0.85, 1.04	0.3	0.91	0.81, 1.01	0.076
Baseline LDH	370	0.98	0.90, 1.06	0.5	0.97	0.89, 1.06	0.5
ECOG 2+	370	1.08	0.70, 1.68	0.7	0.97	0.60, 1.54	0.9
True EMD	370	0.58	0.33, 1.02	0.061	0.53	0.29, 0.94	0.033
Prior LOTs	370	0.92	0.84, 1.00	0.045	0.87	0.79, 0.96	0.006
Prior BCMA	370	1.17	0.77, 1.76	0.5	1.84	1.12, 3.04	0.016
Steroid use
BsAb	355
Tec		—	—		—	—
Elra		0.83	0.48, 1.40	0.5	0.87	0.50, 1.50	0.6
Baseline Hgb	355	0.89	0.79, 1.01	0.074	0.90	0.79, 1.03	0.14
Baseline LDH	355	1.07	0.98, 1.17	0.12	1.05	0.96, 1.15	0.3
ECOG 2+	355	1.71	1.03, 2.83	0.036	1.32	0.78, 2.23	0.3
True EMD	355	1.61	0.85, 2.97	0.13	1.57	0.81, 2.94	0.2
Prior LOTs	355	0.88	0.79, 0.98	0.028	0.94	0.83, 1.06	0.3
Prior BCMA	355	0.47	0.29, 0.78	0.003	0.65	0.36, 1.16	0.15
1 Estimate = Odds ratio or hazard ratio

8 NON-LINEAR ANALYSIS

8.1 OS based on age

data1 <- IPTW.data %>%
  transmute(
    site,
    age_at_bs_ab_initiation,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )


dd <- datadist(data1)
options(datadist='dd')

ref_value = 85

S <- Surv(data1$Days.to.death.or.DLC, data1$Death)
f <- cph(S ~ rcs(age_at_bs_ab_initiation,3) * bs_ab_name , x=TRUE, y=TRUE,surv=TRUE,data=data1, weights = weights)

coxph( S ~ age_at_bs_ab_initiation * bs_ab_name , data=data1) %>% tbl_regression(exponentiate=TRUE)

Characteristic	HR1	95% CI1	p-value
age_at_bs_ab_initiation	0.99	0.97, 1.01	0.3
bs_ab_name
Tec	—	—
Elra	0.05	0.00, 1.01	0.051
age_at_bs_ab_initiation * bs_ab_name
age_at_bs_ab_initiation * Elra	1.05	1.00, 1.09	0.029
1 HR = Hazard Ratio, CI = Confidence Interval

model <- Predict(f, bs_ab_name, age_at_bs_ab_initiation, fun=function(x) exp(x) /exp(predict(f, data.frame(age_at_bs_ab_initiation = ref_value, bs_ab_name = "Tec")) )    )

OS_age <- ggplot(as.data.frame(model),aes(x=age_at_bs_ab_initiation, y=yhat, z = bs_ab_name)) + 
  geom_line(
    aes(col = bs_ab_name)
  ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Elra"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "steelblue" ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Tec"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "coral2" ) + 
  theme_bw() + 
  labs(
    title = "A)",
    x = "Age (years)",
    y = "OS (HR)") +
  scale_color_manual(values = c ("coral2", "steelblue")) + 
  theme(title = element_text(size=16,face="bold"),legend.title=element_blank()) + 
  #coord_cartesian(ylim = c(0.5, 5.5)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray", size = 0.5)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

OS_age

8.2 PFS based on age

data1 <- IPTW.data %>%
  transmute(
    site,
    age_at_bs_ab_initiation,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )


dd <- datadist(data1)
options(datadist='dd')

ref_value = 85

S <- Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death)
f <- cph(S ~ rcs(age_at_bs_ab_initiation,3) * bs_ab_name , x=TRUE, y=TRUE,surv=TRUE,data=data1, weights = weights)

coxph( S ~ age_at_bs_ab_initiation * bs_ab_name , data=data1) %>% tbl_regression(exponentiate=TRUE)

Characteristic	HR1	95% CI1	p-value
age_at_bs_ab_initiation	0.98	0.97, 1.00	0.034
bs_ab_name
Tec	—	—
Elra	0.54	0.06, 4.75	0.6
age_at_bs_ab_initiation * bs_ab_name
age_at_bs_ab_initiation * Elra	1.01	0.98, 1.04	0.4
1 HR = Hazard Ratio, CI = Confidence Interval

model <- Predict(f, bs_ab_name, age_at_bs_ab_initiation, fun=function(x) exp(x) /exp(predict(f, data.frame(age_at_bs_ab_initiation = ref_value, bs_ab_name = "Tec")) )    )

PFS_age <- ggplot(as.data.frame(model),aes(x=age_at_bs_ab_initiation, y=yhat, z = bs_ab_name)) + 
  geom_line(
    aes(col = bs_ab_name)
  ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Elra"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "steelblue" ) + 
  geom_ribbon(data = filter(model, bs_ab_name =="Tec"), aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0, fill = "coral2" ) + 
  theme_bw() + 
  labs(
    title = "B)",
    x = "Age (years)",
    y = "PFS (HR)") +
  scale_color_manual(values = c ("coral2", "steelblue")) + 
  theme(title = element_text(size=16,face="bold"),legend.title=element_blank()) + 
  #coord_cartesian(ylim = c(0.5, 5.5)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray", size = 0.5)

PFS_age

8.3 ***FIGURE 5: OS/PFS and age

combined_age <- (OS_age / PFS_age) +
  plot_layout(guides = "collect", axes = "collect")


ggsave("svg/FIGURE5.svg", combined_age, device = "svg",
       width = 6.5,
       height = 6.5,
       units = c("in"))

ggsave("jpeg/FIGURE5.jpeg", combined_age, device = "jpeg",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE5.tiff", combined_age, device = "tiff",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE5.jpeg")

9 CONTOUR PLOTS

9.1 FIGURE 6: Baseline Hgb and LDH

data1 <- IPTW.data %>%
  filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) %>% ## Exclude patients who received bsAb as bridging
  #filter(bs_ab_name == "tec") %>%
  transmute(
    bs_ab_name,
    weights,
    "CR/sCR" = ifelse(best_response %in% c("CR", "sCR"),1,0),
    "PD" = ifelse(best_response %in% c("PD"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    ECOG_2_4 = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    baseline_hgb,
    baseline_ldh = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    baseline_ferritin = as.numeric(baseline_ferritin_prior_to_bs_ab_initiation),
    ) %>% filter(!is.na(baseline_ldh))


## Tec
data_tec = data1 %>% filter(bs_ab_name == "tec")

dd_tec <- datadist(data_tec)
options(datadist='dd_tec')

S_tec <- Surv(data_tec$Days.to.death.or.relapse.or.DLC/30, data_tec$Relapse.or.death)
f_tec <- cph(S_tec ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_tec, weights=weights)
med_tec <- Quantile(f_tec)

pred_tec  <- data.frame(Predict(f_tec,  baseline_hgb, baseline_ldh,
                                fun = function(x) med_tec(lp = x))) %>%
  mutate(drug = "Teclistamab")

## Elra
data_elra = data1 %>% filter(bs_ab_name == "elra")

dd_elra <- datadist(data_elra)
options(datadist='dd_elra')

S_elra <- Surv(data_elra$Days.to.death.or.relapse.or.DLC/30, data_elra$Relapse.or.death)
f_elra <- cph(S_elra ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_elra, weights=weights)
med_elra <- Quantile(f_elra)

pred_elra <- data.frame(Predict(f_elra, baseline_hgb, baseline_ldh,
                                fun = function(x) med_elra(lp = x))) %>%
  mutate(drug = "Elranatamab")

## Combine both predictions
all_pred <- bind_rows(pred_tec, pred_elra)
z_range <- range(all_pred$yhat, na.rm = TRUE)
brks    <- pretty(z_range, n = 10)   # tweak n if you want more/fewer bands

all_pred <- all_pred %>%
  mutate(
    z_band = cut(
      yhat,
      breaks        = brks,
      include.lowest = TRUE,
      right          = FALSE
    )
  )

z_levels <- levels(all_pred$z_band)

cols <- viridis_pal()(length(z_levels))


contour_plot_tec <- all_pred %>%
  filter(drug == "Teclistamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +   # or geom_tile(); blocky but guarantees alignment
  theme_classic() +
  theme(
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
  )+
  labs(
    title = "A) Teclistamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)",
    fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

contour_plot_tec

contour_plot_elra <- all_pred %>%
  filter(drug == "Elranatamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +
  theme_classic() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
    ) +
  labs(
    title = "B) Elranatamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)"
    #fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

 contour_plot_elra

9.2 ***FIGURE 6: Hgb and LDH contour

combined_contour <- (contour_plot_tec / contour_plot_elra) +
  plot_layout(guides = "collect", axes = "collect")


ggsave("svg/FIGURE6.svg", combined_contour, device = "svg",
       width = 6.5,
       height = 6.5,
       units = c("in"))

ggsave("jpeg/FIGURE6.jpeg", combined_contour, device = "jpeg",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE6.tiff", combined_contour, device = "tiff",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE6.jpeg")